function pliMCmain(inDir,outDir,outAscii)

%inDir='C:\Users\matthew\Desktop\Projects\Mental-Calculation\data\';
%outDir='C:\Users\matthew\Desktop\Projects\Batch\analysis\mc\';
%outAscii='C:\Users\matthew\Desktop\Projects\Batch\analysis\mc\ascii\';

pliXBand(inDir,outDir,outAscii);
%pliStat(outDir,outDir);
function pliXBand(inDir,outDir,outAscii)

sampleOI=4096;
nTrials=1;
bands={[1 4],[4 8],[8 13],[13 30],[30 49]};
bandLabels={'delta','theta','alpha','beta','gamma'};
nBands=length(bandLabels);



filter=fullfile(inDir,'*comp*.ds');
files=dir(filter);
nSubjs=length(files);


for i=1:nSubjs
    subjName=files(i).name;
    trials=trial(fullfile(inDir,subjName))
    trials=chunkOI(trials,sampleOI,nTrials)
    
    for j=1:nBands   
        
        bf_trials=selectBand(trials,bands{j});
        
        asciiFile=strcat(bandLabels{j},'_');
        subjName=strtok(subjName,'.');
        subjName=strtok(subjName,'0');
        asciiFile=strcat(asciiFile,subjName);
        
        %write2ascii(fullfile(outAscii,strcat(asciiFile,'.txt')),bf_trials,nTrials,precision);
        
        [pliXT ,spliXT]=pliXTrial(bf_trials);
        
        outfile=strcat('pli_',bandLabels{j});
        outfile=strcat(outfile,'_');
        
        outfile=strcat(outfile,subjName);
        
        save(fullfile(outDir,outfile),'pliXT','spliXT');
    end    
end


function [trials]=trial(filename)

% Band filtering
cfg= [];
% Set data
cfg.dataset= filename;
% Remove bad channels, set the trigger, set pre and post-stimulus interval
cfg.channel = {'MEG','-MLT22','-MLO42','-MLO43','-MRT22','-MRO42','-MRO43','-MZ'};
%Call definetrial to compute trials



cfg.headerformat = 'ctf_old';
cfg.dataformat= 'ctf_old';

cfg.continuous='yes';
trials=ft_preprocessing(cfg);

function [trials]=chunkOI(trials,sampleOI,nTrials)
%trials.trial{1}=trials.trial{1}(:,sampleOI);
%trials.time{1}=trials.time{1}(1,sampleOI);
cfg=[];
cfg.length=(sampleOI*1/trials.fsample)%%%occhio!!!  
cfg.overlap=0;
trials=ft_redefinetrial(cfg,trials);
trials.trial=trials.trial(1:nTrials);
trials.time=trials.time(1:nTrials);
for i=1:length(trials.trial)
    trials.trial{i}=trials.trial{i};
end



function [bf_trials]=selectBand(trials,band)

bf_trials=trials;
nTrials=length(trials.trial);
cfg = [];
cfg.bpfilter='yes';
cfg.bpfreq=band;
bf_trials= ft_preprocessing(cfg,trials);
% for i=1:nTrials
%     
%     bf_trials.trial{i}=eegfilt(trials.trial{i},trials.fsample,band(1),band(2));
% 
% end
function [pli3d,spli3d]=pliXTrial(bf_trials)

nTrials=length(bf_trials.trial);
nChs=size(bf_trials.trial{1},1);
pli3d=zeros(nTrials,nChs,nChs);
spli3d=zeros(nTrials,nChs,nChs);
for i=1:nTrials
    [pli3d(i,:,:),zscorematrix,spli3d(i,:,:)]=pli2(bf_trials.trial{i}');  
end

% function PLI=pli(a)
% 
% % a is a filtered multichannel signal (time x channels) hilbert(a)
% % calculates analytic signal (complex valued) of each column of a. Phase
% % Lag Index between channel i and j averaged over time bins is stored in
% % PLI(i,j) number of channels
% N=size(a,2);
% PLI(1:N,1:N)=0;
% complex_a=hilbert(a);
% for i=1:N-1
%    for j=i+1:N
%         PLI(i,j)=abs(mean(sign(imag(complex_a(:,i)./complex_a(:,j)))));
%    end
% end
% PLI=PLI+PLI';

%change values of areas.. they are wrong
function pliStat(inDir,outDir)
mlc=1:15;
mlf=16:31;
mlo=32:39;
mlp=40:48;
mlt=49:68;

mrc=69:83;
mrf=84:99;
mro=100:107;
mrp=108:116;
mrt=117:136;

labels={'mlc','mlf','mlo','mlp','mlt','mrc','mrf','mro','mrp','mrt'};
areas={mlc;mlf;mlo;mlp;mlt;mrc;mrf;mro;mrp;mrt};
bandlabels={'theta','alpha','beta','gamma'};
alphalevel=0.05;
for i=1:length(bandlabels)
    
    filter=strcat('pli_',bandlabels{i});
    filter=strcat(filter,'_sp*')
    subjs=dir(fullfile(inDir,filter));
    for j=1:length(subjs)
        
        subjName=subjs(j).name;
        load(fullfile(inDir,subjName));
        one_pliXT=pliXT;
        
        [r,subjName]=strtok(subjName,'-');
        
        subjName=strcat(strcat(strcat('pli_',bandlabels{i}),'_comp'),subjName);
        load(fullfile(inDir,subjName));
        
        two_pliXT=pliXT;

        statsingle=statOverTr(one_pliXT,two_pliXT,areas,alphalevel);
        [r,subjName]=strtok(subjName,'-');
        outfile=strcat('stat_',bandlabels{i});
        outfile=strcat(outfile,subjName);
        save(fullfile(outDir,outfile),'statsingle');
    end
    
    statOverSubjs(inDir,outDir,bandlabels{i},alphalevel);
    
end


function [res]=statOverTr(one_pliXT,two_pliXT,areas,alphalevel)

nTrial=min(size(one_pliXT,1),size(two_pliXT,1));
one_trialXintra_area=zeros(nTrial,length(areas));
two_trialXintra_area=zeros(nTrial,length(areas));

interAreaLen=((length(areas)*length(areas))-length(areas))/2;

one_trialXinter_area=zeros(nTrial,interAreaLen);
two_trialXinter_area=zeros(nTrial,interAreaLen);


s=size(one_pliXT);
indx=ones(s(2),s(3));
indx=indx-tril(indx);
indx=logical(indx);

one_avg_pliXtrial=zeros(nTrial,1);
two_avg_pliXtrial=zeros(nTrial,1);


for i=1:nTrial
   
   tone=squeeze(one_pliXT(i,:,:));
   ttwo= squeeze(two_pliXT(i,:,:));
    
   one_trialXintra_area(i,:)=shortPLI(tone,areas);
   two_trialXintra_area(i,:)=shortPLI(ttwo,areas);
   
   one_trialXinter_area(i,:)=longPLI(tone,areas);
   two_trialXinter_area(i,:)=longPLI(ttwo,areas);
   
   one_avg_pliXtrial(i)=mean(tone(indx));
   two_avg_pliXtrial(i)=mean(ttwo(indx));

end

res.intra={one_trialXintra_area,two_trialXintra_area};
res.inter={one_trialXinter_area,two_trialXinter_area};
res.avg={one_avg_pliXtrial,two_avg_pliXtrial};

[res.h_intra,res.p_intra]=ttest(one_trialXintra_area,two_trialXintra_area,alphalevel);
[res.h_inter,res.p_inter]=ttest(one_trialXinter_area,two_trialXinter_area,alphalevel);
[res.h_avgpli,res.p_avgpli]=ttest(one_avg_pliXtrial,two_avg_pliXtrial,alphalevel);
function statOverSubjs(inDir,outDir,bandName,alphalevel)

filter=strcat('stat_*',bandName);
filter=strcat(filter,'*');

files=dir(fullfile(inDir,filter));
nSubjs=length(files);

load(fullfile(inDir,files(1).name),'statsingle');

nAreas=size(statsingle.intra{1},2)
nPairAreas=size(statsingle.inter{1},2);

one_SubjXintraA=zeros(nSubjs,nAreas);
two_SubjXintraA=zeros(nSubjs,nAreas);
one_SubjXinterA=zeros(nSubjs,nPairAreas);
two_SubjXinterA=zeros(nSubjs,nPairAreas);
one_SubjAVG=zeros(nSubjs,1);
two_SubjAVG=zeros(nSubjs,1);

for i=1:nSubjs
    load(fullfile(inDir,files(i).name),'statsingle');

    one_SubjXintraA(i,:)=mean(statsingle.intra{1},1);
    two_SubjXintraA(i,:)=mean(statsingle.intra{2},1);
    one_SubjXinterA(i,:)=mean(statsingle.inter{1},1);
    two_SubjXinterA(i,:)=mean(statsingle.inter{2},1);
    one_SubjAVG(i)=mean(statsingle.avg{1},1);
    two_SubjAVG(i)=mean(statsingle.avg{2},1);
    
    [r,subjName]=strtok(files(i).name,'-');
    subjName=subjName(2:end);
    subjName=strtok(subjName,'.');
    multi.subjNames{i}=subjName;  
end

multi.intra={one_SubjXintraA,two_SubjXintraA};
multi.inter={one_SubjXinterA,two_SubjXinterA};
multi.avg={one_SubjAVG,two_SubjAVG};

[multi.h_intra,multi.p_intra]=ttest(one_SubjXintraA,two_SubjXintraA,alphalevel);
[multi.h_inter,multi.p_inter]=ttest(one_SubjXinterA,two_SubjXinterA,alphalevel);
[multi.h_avgPLI,multi.p_avgPLI]=ttest(one_SubjAVG,two_SubjAVG,alphalevel);

outfile=fullfile(outDir,strcat('statMultiSubjs_',bandName));

save(outfile,'multi');
function [v]=shortPLI(pli,areas)

v=zeros(1,length(areas));
for i=1:length(areas)
    m=pli(areas{i},areas{i}); 
    indx=ones(size(v));
    indx=indx-tril(indx);
    indx=logical(indx);    
    v(i)=mean(m(indx));
end
function [interArea]=longPLI(pli,areas)

v=zeros(length(areas),length(areas));
for i=1:length(areas)-1
    for j=i+1:length(areas)
        m=pli(areas{i},areas{j});
        s=size(m);
        v(i,j)=mean(reshape(m,s(1)*s(2),1));
    end
end

indx=ones(size(v));
indx=indx-tril(indx);
indx=logical(indx);
%indx=logical(indx);%%warning if phase is 0
interArea=v(indx)';


function [PLI,zscoreMatrix,sPLI]=pli2(m)

% a is a filtered multichannel signal (time x channels)
% hilbert(a) calculates analytic signal (complex valued) of each
% column of a. Phase Lag Index between channel i and j averaged over
% time bins is stored in PLI(i,j)
% number of channels

nChs=size(m,2);
timeStop=size(m,1);

PLI=zeros(nChs,nChs);
complex_m=hilbert(m);

phi=atan2(imag(complex_m),real(complex_m));


zscoreMatrix=zeros(nChs,nChs);

for i=1:nChs
    
    for j=1:nChs
        count=0;
        count2=0;

        for t=1:timeStop
        
            phaseDiff=phi(t,i)-phi(t,j);
            
            if(sin(phaseDiff)>0)
                count=count+1;
            end
            if(abs(sin(phaseDiff))>0)
                count2=count2+1;
            end
            
        end
        
        if(count2>0)
            pli=2*abs(.5-1*count/count2);
            %for FDR
            zscoreMatrix(i,j)=((count/count2)-0.5)/sqrt((0.5*0.5)/count2);
        else
            pli=0;
        end
        
        PLI(i,j)=pli;
    end
end
sPLI=fdr(zscoreMatrix,PLI,0.001);
function [pliM]=fdr(zscoreMatrix,pliM,alpha)

lindx=find(triu(ones(size(pliM)),1));
[row,col]=ind2sub(size(pliM),lindx);
z=abs(zscoreMatrix(lindx));
k=ones(1,length((z)));%zeros or ones?? stam one
N=length(z);

for i=1:length(z)
    
    k(i)=k(i)+sum(z(i)<z);%to be equal of stam code
    
    if(z2p(z(i))>(k(i)*alpha)/N)
        pliM(row(i),col(i))=0;
        pliM(col(i),row(i))=0;
    end        
        
end
function [p]=z2p(zscore)
        p=1;
       
        if(zscore > 0.67 )
            p=0.5;
        end
        if(zscore > 0.84)
            p=0.4;
        end
        if(zscore > 1.04 )
            p=0.3;
        end
        if(zscore > 1.28 )
            p=0.2;
        end
        if(zscore > 1.64)
            p=0.1;
        end
        if(zscore > 1.96 )
            p=0.05;
        end
        if(zscore > 2.33 )
            p=0.02;
        end
        if(zscore > 2.58 )
            p=0.01;
        end
        if(zscore > 2.81 )
            p= 0.005;
        end
        if(zscore > 3.09 )
            p=0.002;
        end
        if(zscore > 3.29 )
            p=0.001;
        end
        if(zscore > 3.89 )
            p=0.0001 ;
        end


